library(tidyverse)
library(skimr)
library(rethinking)

Notes

Introduction

Linear regression describes a measurements mean and variance as the addition of other measurements. It assumes the errors in the primary measurement are of the gaussian distribution.

Normality

1000 individuals conduct a random walk on either side of the football field starting at the 50 yard line (value = 0). This signifies the normal distribution.

INDIVIDUALS <- 1000
TRIALS <- 100

purrr::map(1:INDIVIDUALS, ~ runif(TRIALS, -1, 1)) %>% 
  purrr::map(cumsum) %>% 
  reduce(rbind) %>% 
  as_tibble %>% 
  `colnames<-`(paste0(1:TRIALS, "t")) %>% 
  mutate(individual = row_number()) %>% 
  reshape2::melt("individual") %>% 
  arrange(individual, variable) %>% 
  ggplot(aes(x=variable, y=value, group=individual)) + 
    geom_line(alpha = I(1/20), size=1) +
    geom_boxplot(aes(x=variable, y = value), inherit.aes = FALSE, colour = "yellow", alpha= I(1/40))

The above example shows us that normality is additive.

# Generate 12 numbers between 1 and 1.1 then add them in to one sum
# Do this 1000 times and collect the numbers in a vector
growth <- replicate(1000, sum(runif(12,0,.1)))
# Plot the distribution of thes enumbers. Note that it is normal.
qplot(growth)

We can also show that normality is multiplicative.

# Generate 12 numbers between 1 and 1.1 then multiply them in to one product
# Do this 1000 times and collect the numbers in a vector
growth <- replicate(1000, prod(1 + runif(12,0,.1)))
#plot the distribution of these numbers. Note that it is normal.
qplot(INDIVIDUALS)

Things start to skew at very large deviates, though this can be corrected by log scaling.

big <- replicate(10000, prod(1 + runif(12, 0, .5)))
small <- replicate(10000, prod(1 + runif(12, 0, 0.01)))

data_frame(big = big, 
           small = small, 
           log_big = log(big)) %>% 
  gather %>% 
  ggplot(aes(value)) + 
    geom_histogram() + 
    facet_wrap("key", scales = "free")

A language for describing models

Elements:

An outcome variable: The set of measurements we hope to predict or understand
Likeliehood distribution: For each of these outcome variables, we define the plausibility of that observations.
Predictor variables: Set of other measurements that we hope to use to predict or understand the outcome.

Approach:

Relate the exact shape of the likelihood distribution to the predictor variables via parameters
Choose priors for all the parameters in the model.

Example: A gaussian model of height

We’ll build a regression to predict a human’s hight, which we will model as a Gaussian distribution with two parameters, mean and standard deviation. There are an infinite number of possible Gaussian distributions. We want our Bayesian machine to consider every possible distribution, each defined by a combination of meand and sd, and rank them by posterior plausibility given the data.

data(Howell1)
d <- Howell1
d2 <- d %>% 
  filter(age >= 18)
d2 %>% 
  skim

Skim summary statistics
 n obs: 352 
 n variables: 4 

── Variable type:integer ─────────────────────────────────────────────────────────
 variable missing complete   n mean  sd p0 p25 p50 p75 p100     hist
     male       0      352 352 0.47 0.5  0   0   0   1    1 ▇▁▁▁▁▁▁▇

── Variable type:numeric ─────────────────────────────────────────────────────────
 variable missing complete   n   mean    sd     p0    p25    p50    p75   p100
      age       0      352 352  41.14 15.97  18     28     39     51     88   
   height       0      352 352 154.6   7.74 136.53 148.59 154.31 160.66 179.07
   weight       0      352 352  44.99  6.46  31.07  40.26  44.79  49.29  62.99
     hist
 ▇▇▇▅▂▃▁▁
 ▁▅▇▆▆▃▁▁
 ▂▆▇▇▇▃▂▁

Lets specifiy our model like so. The first line is the likelihood of our outcome, height. The second two lines are priors for the parameters in our likelihood. \[ h_i \sim Normal(\mu, \sigma) \\ \mu \sim Normal(178, 20) \\ \sigma \sim Uniform(0, 50) \]

As height is a physical phenomenon, we can use our intuition to assign a prior.

data_frame(x=100:250, 
           y=purrr::map_dbl(100:250, ~ dnorm(.x, 178, 20))) %>% 
  ggplot(aes(x,y)) + geom_path()

We give the standard deviation a truly flat prior. We really just want to make sure its positive.

data_frame(x=-10:60, 
           y=purrr::map_dbl(-10:60, ~ dunif(.x, 0, 50))) %>% 
  ggplot(aes(x,y)) + geom_path()

Specifying a prior for the parameters of our outcome implicitly give us a prior for our outcome.

sample_mu <- rnorm(1e4, 178, 20)
sample_sigma <- runif(1e4, 0, 50)
prior_h <- rnorm(1e4, sample_mu, sample_sigma)
qplot(prior_h)

Lets us grid approximation to generate the posterior. To do so, we enumerate a list of candidate values for our parameters and store them in post. We then calculate the log likelihood for each of those parameters given all the heights we have seen and store them in the column LL. We use log transformation to avoid rounding errors. Then we multiply by the priors (log addition is the same as multiplying)

# Generate candidate parameters
mu_list <- seq(from=140, to=160, length.out=200)
sigma_list <- seq(from=4, to=9, length.out=200)
post <- expand.grid(mu=mu_list, sigma=sigma_list)

# Calculate the likelihood of every possible
# combination candidate parameters
post$LL <- sapply(1:nrow(post), function(i) {
  sum(dnorm(
    d2$height,
    mean=post$mu[i],
    sd=post$sigma[i],
    log=TRUE))
  })

# Multiply priors
post$prod <- post$LL + dnorm(post$mu, 178, 20, TRUE) +
  dunif(post$sigma, 0, 50, TRUE)

# Calculated scaled posterior
post$prob <- exp(post$prod - max(post$prod))

# Contour plot of probabilities
ggplot(post, aes(x=mu, y=sigma, z=prob)) + 
  geom_contour(aes(colour=stat(level))) + 
  coord_cartesian(xlim=c(140, 160), ylim=c(4, 9))

Now lets sample from the posterior.

# Sample candidate values by posterior
sample_rows <- sample(1:nrow(post), size=1e4, replace=TRUE, prob=post$prob)

# Plot their 2d histogram
data_frame(
  sample_mu = post$mu[sample_rows],
  sample_sigma = post$sigma[sample_rows]
) %>% 
  ggplot(aes(sample_mu, sample_sigma)) +
    geom_bin2d(bins=30) + 
    coord_cartesian(xlim=c(140, 160), ylim=c(4, 9))

Lets also examine the margins.

require(gridExtra)

Loading required package: gridExtra

Attaching package: ‘gridExtra’

The following object is masked from ‘package:dplyr’:

    combine

p1 <- post$mu[sample_rows] %>% 
  qplot + ggtitle("mu")
p2 <- post$sigma[sample_rows] %>% 
  qplot + ggtitle("sigma")

grid.arrange(p1, p2, ncol=2)

map from the rethinking package allows us to solve for the posterior using quadratic approximation instead of grid approximations. It shows us gaussian approximations for each parameter’s marginal distribution.

flist <- alist(
  height ~ dnorm(mu, sigma),
  mu ~ dnorm(178, 20),
  sigma ~ dunif(0, 50)
)
m4.1 <- rethinking::map(flist, data=d2)
precis(m4.1)

Variance covariance matrix tells us how each parameter relates to every other parameter in the posterior distribution.

vcov(m4.1)

                mu        sigma
mu    0.1697397726 0.0002182674
sigma 0.0002182674 0.0849060205

# This can be factored in to a vector of variances
diag(vcov(m4.1))

        mu      sigma 
0.16973977 0.08490602

# Or a correlation matrix
cov2cor(vcov(m4.1))

               mu       sigma
mu    1.000000000 0.001818142
sigma 0.001818142 1.000000000

Lets draw from the posterior by sampling vectors of values from a multi-dimensional gaussian distribution.

post <- extract.samples(m4.1, n=1e4)
precis(post)

Adding a predictor variable, weight

It looks like weight will be a good predictor for height in adults.

ggplot(d2, aes(x=weight, y=height)) + 
  geom_point(shape=1) +
  geom_smooth(method=lm)

The goal here is to make the parameter for the mean of a Gaussian distribution, \(\mu\), into a linear function of the predictor variable. This encodes the assumption that the predictor variable has a perfectly constant and additive relationship to the mean of the outcome.

m4.3 <- rethinking::map(
  alist(
    height ~ dnorm(mu, sigma),
    mu <- a + b*weight,
    a ~ dnorm(156, 100),
    b ~ dnorm(0, 10),
    sigma ~ dunif(0, 50)),
  data = d2)
# Display table of estimates
precis(m4.3, corr=TRUE)

The table above tells us that a person 1 kg heavier is expected to be .9 cm taller. It also tells us that a person with 0 weight is 114 cm tall, which is of course nonsense; This is why it is usually important to have very weak priors for intercepts in many cases. Finally, it explains that 95% of plausible heights lie within 10 cm of the mean.

Centering is the procedure of subtracting the mean of a variable from each value.

d2$weight.c <- d2$weight - mean(d2$weight)

m4.4 <- rethinking::map(
  alist(
    height ~ dnorm(mu, sigma),
    mu <- a + b*weight.c,
    a ~ dnorm(178, 100),
    b ~ dnorm(0, 10),
    sigma ~ dunif(0, 50)),
  data = d2)
# Display table of estimates
precis(m4.4, corr=TRUE)

By centering our predictor variable, the interpretation of the intercept has now become the expected value of the outcome when the predictor is at its average value. This makes interpreting the intercept a lot easier.

Visualizing our estimates

We can plot our maximum a posteriori fit here

p <- ggplot(d2, aes(x=weight, y=height)) +
  geom_point(shape=1, size=2) +
  geom_abline(aes(intercept = coef(m4.3)["a"],
                  slope = coef(m4.3)["b"]))
p

This line only represents one piece of the posterior. We could instead sample values of \(\alpha\) and \(\beta\) from the posterior to generate many lines to build bounds of uncertainty. We’ll do this with a much smaller dataset to magnify effects.

N <- 10
dN <- d2[ 1:N , ]
mN <- rethinking::map(
  alist(
    height ~ dnorm( mu , sigma ) ,
    mu <- a + b*weight ,
    a ~ dnorm( 178 , 100 ) ,
    b ~ dnorm( 0 , 10 ) ,
    sigma ~ dunif( 0 , 50 )
) , data=dN )

post <- extract.samples(mN, n=20)
head(post)

p <- ggplot(dN, aes(weight, height)) +
  geom_point(shape=1)

for (i in 1:20) {
  p <- p + geom_abline(intercept = post$a[i], slope = post$b[i], alpha = I(3/10))
}

p

To do this more generally:

Sample parameter values from the posterior.
Use the parameter values to create predictions of the parameters of interest, including the final outcome.
Use summary functions like mean, HDOP, PI to find averages and lower and upper bounds of our uncertainty.

# Draw parameter values from the posterior distribution
posterior_samples <- extract.samples(m4.3)

# Create out x-axis of interest
weights_axis <- 25:70

# Function to calculate fit values for \mu
mu_link <- function(weight) post$a + post$b*weight
mu <- purrr::map(weights_axis, mu_link)

# Function to calculate fit heights
# This incorporates variability as well
simulate <- function(weight, parameter_samples) {
  rnorm(n = nrow(parameter_samples),
        mean = parameter_samples$a + parameter_samples$b * weight,
        sd = parameter_samples$sigma)}
simulated_heights <- purrr::map(weights_axis, ~ simulate(.x, posterior_samples))

# Combine all this information to plot
data_frame(
  # X axis
  weight = weights_axis,
  # Estimates for mean height
  mu_mean = map_dbl(mu, mean),
  mu_lower = map_dbl(mu, ~HPDI(.x, .89)[1]),
  mu_upper = map_dbl(mu, ~HPDI(.x, .89)[2]),
  # Estimates for height
  height_mean = map_dbl(simulated_heights, mean),
  height_low = map_dbl(simulated_heights, ~HPDI(.x, .89)[1]),
  height_high = map_dbl(simulated_heights, ~HPDI(.x, .89)[2])) %>% 
  ggplot(data=.) +
  geom_point(inherit.aes = FALSE, data = d2,
             aes(x=weight, y=height), shape=1) +
  geom_line(aes(x=weight, y=mu_mean)) +
  geom_ribbon(aes(x = weight, ymin=mu_lower, ymax=mu_upper), alpha=.5) +
  geom_ribbon(aes(x = weight, ymin = height_low, ymax = height_high), alpha=.3) +
  ggtitle("The effect of weight on height", subtitle = "Using 89% highest posterior density intervals") +
  coord_cartesian(xlim=c(30, 60),ylim=c(135, 180))

Polynomial Regression

Standardizing our predictors means scaling but also dividing by the standard deviation. This changes our interpretation of continuous variables to mean the unit change in outcome from a cahnge of one standard deviation.

# Standardize continuous variable
standardize <- function(vec) {
  centered <- vec - mean(vec)
  centered / sd(vec)
}
d$weight.s <- standardize(d$weight)

# Define another polynomial term
d$weight.s2 <- d$weight.s^2

m4.5 <- map(
  alist(
    # Outcome
    height ~ dnorm(mu, sigma),
    # Regression specification
    mu <- a + b1*weight.s + b2*weight.s2,
    a ~ dnorm(178 , 100),
    # These priors are very weak!
    b1 ~ dnorm(0 , 10),
    b2 ~ dnorm(0 , 10),
    sigma ~ dunif(0 , 50)),
  data=d )
precis( m4.5 )

On transformations

Centering and scaling a continuous variable increases interpretability. Demeaning a predictor variable subtracts every value by it’s mean. By demeaning a predictor variables, the intercept represents the value of the outcome when all variables are at their mean. One standardizes a variable by dividing every value by the variables standard deviation. Standardizing causes the interpretation of coefficients to be in terms of unit changes in SD of the variable.

Log transforming the outcome variable makes y our coefficients represent percent increases in the outcome. Log transforming a predictor variable describes how a single percentage point increase in the predictor affects y. Log transforming both the predictor and outcome describes how percentage changes in one create percentage changes in the other.

Homework

Medium

4M1

For the model definition below, simulate observed heights from the prior.

num_observations <- 1e4

prior_mean <- rnorm(num_observations, 0, 10)
prior_sd <- runif(num_observations, 0, 10)

prior_heights <- rnorm(num_observations, prior_mean, prior_sd)

qplot(prior_heights)

4M2

model <- alist(
  y ~ dnorm(mu, sigma),
  mu ~ dnorm(0, 10),
  sigma ~ dunif(0, 10)
)

Hard

4H1

data("Howell1")
d <- Howell1
d$weight <- standardize(d$weight)
d

model <- map(
  alist(
    height ~ dnorm(mu, sigma),
    mu <- a + b * weight,
    a ~ dnorm(mean = 0, sd = 100),
    b ~ dnorm(mean = 0, sd = 10),
    sigma ~ dunif(min = 0, max = 64)
  ),
  data = d
)
model


Maximum a posteriori (MAP) model fit

Formula:
height ~ dnorm(mu, sigma)
mu <- a + b * weight
a ~ dnorm(mean = 0, sd = 100)
b ~ dnorm(mean = 0, sd = 10)
sigma ~ dunif(min = 0, max = 64)

MAP values:
         a          b      sigma 
138.261398  25.927300   9.345958 

Log-likelihood: -1987.71

# Standardize our input
new_weights <- standardize(c(46.95, 43.72, 64.78, 32.59, 54.63))

# Draw parameter values from the posterior distribution
posterior_samples <- extract.samples(model)

# Function to calculate fit heights
simulate_heights <- function(weight, parameter_samples) {
  rnorm(
    n = nrow(parameter_samples),
    mean = parameter_samples$a + parameter_samples$b * weight,
    sd = parameter_samples$sigma)}

# Generate a prediction for each new weight
# Use every parameter set in the posterior to do so
simulated_heights <- purrr::map(new_weights, ~ simulate(.x, posterior_samples))

data_frame(
  individual = 1:5,
  weight = c(46.95, 43.72, 64.78, 32.59, 54.63),
  expected_height = map_dbl(simulated_heights, mean),
  lower_pi = map_dbl(simulated_heights, ~PI(.x)[1]),
  upper_pi = map_dbl(simulated_heights, ~PI(.x)[2])
)

4H2

data("Howell1")
d <- Howell1
d %>% 
  filter(age < 18) %>% 
  ggplot(aes(x = weight, y = height)) +
  geom_point() + geom_smooth(method=lm) +
  geom_smooth(colour = "orange")

d %>% 
  filter(age < 18) %>% 
  map(alist(
    height ~ dnorm(mu, sigma),
    mu <- a + b * weight,
    a ~ dnorm(mean = 100, sd = 100),
    b ~ dnorm(mean = 0, sd = 10),
    sigma ~ dunif(min = 0, max = 50)
  ),
  data = .) ->
model

precis(model)

A one unit increase in weight increases height by 2.72, therefore a 10 unit increase will move height by 27.2.

# Sample the poseterior
posterior_samples <- extract.samples(model)

# Create x-axis
new_weights <- seq(from = 4, to = 45, length.out = 100)

# Calculate the mean
fn_mean <- function(weight) {posterior_samples$a + posterior_samples$b * weight}
fit_means <- purrr::map(new_weights, fn_mean)

# Calculate heights
fn_height <- function(weight) {
  rnorm(
    n = nrow(posterior_samples),
    mean = posterior_samples$a + posterior_samples$b * weight,
    sd = posterior_samples$sigma)
}
fit_heights <- purrr::map(new_weights, fn_height)

data_frame(
  weights = new_weights,
  means = map_dbl(fit_means, mean),
  means_low = map_dbl(fit_means, ~ HPDI(.x, .89)[1]),
  means_high = map_dbl(fit_means, ~ HPDI(.x, .89)[2]),
  heights = map_dbl(fit_heights, mean),
  heights_low = map_dbl(fit_heights, ~ HPDI(.x, .89)[1]),
  heights_high = map_dbl(fit_heights, ~ HPDI(.x, .89)[2]),
) %>% 
ggplot(aes(x = weights)) +
  geom_line(aes(y = means)) +
  geom_ribbon(aes(ymin = means_low, ymax = means_high), alpha = .4) +
  geom_ribbon(aes(ymin = heights_low, ymax = heights_high), alpha = .2) +
  geom_point(inherit.aes = FALSE, data = filter(d, age < 18),
             aes(x=weight, y = height), shape = 1) +
  ylab("height") + xlab("weight") + theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

BRMS

data(Howell1)
d <- Howell1
d2 <- d %>% 
  filter(age >= 18)
d2 %>% 
  skim

Skim summary statistics
 n obs: 352 
 n variables: 4 

── Variable type:integer ─────────────────────────────────────────────────────────
 variable missing complete   n mean  sd p0 p25 p50 p75 p100     hist
     male       0      352 352 0.47 0.5  0   0   0   1    1 ▇▁▁▁▁▁▁▇

── Variable type:numeric ─────────────────────────────────────────────────────────
 variable missing complete   n   mean    sd     p0    p25    p50    p75   p100
      age       0      352 352  41.14 15.97  18     28     39     51     88   
   height       0      352 352 154.6   7.74 136.53 148.59 154.31 160.66 179.07
   weight       0      352 352  44.99  6.46  31.07  40.26  44.79  49.29  62.99
     hist
 ▇▇▇▅▂▃▁▁
 ▁▅▇▆▆▃▁▁
 ▂▆▇▇▇▃▂▁

library(brms)

b4.1 <-
  brm(data = d2, family = gaussian,
      height ~ 1,
      prior = c(prior(normal(178, 20), class = Intercept),
                prior(cauchy(0, 1), class = sigma)),
      iter = 31000, warmup = 30000, chains = 4, cores = 4)

Compiling the C++ model
Start sampling

starting worker pid=44777 on localhost:11970 at 14:01:32.770
starting worker pid=44787 on localhost:11970 at 14:01:33.005
starting worker pid=44797 on localhost:11970 at 14:01:33.233
starting worker pid=44807 on localhost:11970 at 14:01:33.465

SAMPLING FOR MODEL '1190ce989ffb6811cbd6fe4e5d278ba8' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 3.3e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.33 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:     1 / 31000 [  0%]  (Warmup)

SAMPLING FOR MODEL '1190ce989ffb6811cbd6fe4e5d278ba8' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 2e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.2 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:     1 / 31000 [  0%]  (Warmup)

SAMPLING FOR MODEL '1190ce989ffb6811cbd6fe4e5d278ba8' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 2.8e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.28 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:     1 / 31000 [  0%]  (Warmup)
Chain 1: Iteration:  3100 / 31000 [ 10%]  (Warmup)
Chain 2: Iteration:  3100 / 31000 [ 10%]  (Warmup)

SAMPLING FOR MODEL '1190ce989ffb6811cbd6fe4e5d278ba8' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 3.6e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.36 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:     1 / 31000 [  0%]  (Warmup)
Chain 1: Iteration:  6200 / 31000 [ 20%]  (Warmup)
Chain 3: Iteration:  3100 / 31000 [ 10%]  (Warmup)
Chain 2: Iteration:  6200 / 31000 [ 20%]  (Warmup)
Chain 4: Iteration:  3100 / 31000 [ 10%]  (Warmup)
Chain 1: Iteration:  9300 / 31000 [ 30%]  (Warmup)
Chain 3: Iteration:  6200 / 31000 [ 20%]  (Warmup)
Chain 2: Iteration:  9300 / 31000 [ 30%]  (Warmup)
Chain 4: Iteration:  6200 / 31000 [ 20%]  (Warmup)
Chain 1: Iteration: 12400 / 31000 [ 40%]  (Warmup)
Chain 3: Iteration:  9300 / 31000 [ 30%]  (Warmup)
Chain 2: Iteration: 12400 / 31000 [ 40%]  (Warmup)
Chain 4: Iteration:  9300 / 31000 [ 30%]  (Warmup)
Chain 3: Iteration: 12400 / 31000 [ 40%]  (Warmup)
Chain 1: Iteration: 15500 / 31000 [ 50%]  (Warmup)
Chain 4: Iteration: 12400 / 31000 [ 40%]  (Warmup)
Chain 2: Iteration: 15500 / 31000 [ 50%]  (Warmup)
Chain 3: Iteration: 15500 / 31000 [ 50%]  (Warmup)
Chain 1: Iteration: 18600 / 31000 [ 60%]  (Warmup)
Chain 4: Iteration: 15500 / 31000 [ 50%]  (Warmup)
Chain 3: Iteration: 18600 / 31000 [ 60%]  (Warmup)
Chain 2: Iteration: 18600 / 31000 [ 60%]  (Warmup)
Chain 1: Iteration: 21700 / 31000 [ 70%]  (Warmup)
Chain 4: Iteration: 18600 / 31000 [ 60%]  (Warmup)
Chain 3: Iteration: 21700 / 31000 [ 70%]  (Warmup)
Chain 2: Iteration: 21700 / 31000 [ 70%]  (Warmup)
Chain 1: Iteration: 24800 / 31000 [ 80%]  (Warmup)
Chain 4: Iteration: 21700 / 31000 [ 70%]  (Warmup)
Chain 3: Iteration: 24800 / 31000 [ 80%]  (Warmup)
Chain 2: Iteration: 24800 / 31000 [ 80%]  (Warmup)
Chain 1: Iteration: 27900 / 31000 [ 90%]  (Warmup)
Chain 4: Iteration: 24800 / 31000 [ 80%]  (Warmup)
Chain 3: Iteration: 27900 / 31000 [ 90%]  (Warmup)
Chain 1: Iteration: 30001 / 31000 [ 96%]  (Sampling)
Chain 2: Iteration: 27900 / 31000 [ 90%]  (Warmup)
Chain 3: Iteration: 30001 / 31000 [ 96%]  (Sampling)
Chain 1: Iteration: 31000 / 31000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 1.41519 seconds (Warm-up)
Chain 1:                0.058162 seconds (Sampling)
Chain 1:                1.47335 seconds (Total)
Chain 1: 
Chain 4: Iteration: 27900 / 31000 [ 90%]  (Warmup)
Chain 3: Iteration: 31000 / 31000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 1.25799 seconds (Warm-up)
Chain 3:                0.053061 seconds (Sampling)
Chain 3:                1.31106 seconds (Total)
Chain 3: 
Chain 2: Iteration: 30001 / 31000 [ 96%]  (Sampling)
Chain 4: Iteration: 30001 / 31000 [ 96%]  (Sampling)
Chain 2: Iteration: 31000 / 31000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 1.39502 seconds (Warm-up)
Chain 2:                0.049075 seconds (Sampling)
Chain 2:                1.4441 seconds (Total)
Chain 2: 
Chain 4: Iteration: 31000 / 31000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 1.2608 seconds (Warm-up)
Chain 4:                0.040366 seconds (Sampling)
Chain 4:                1.30116 seconds (Total)
Chain 4:

plot(b4.1)

Extract variance covariance matrix.

post <- posterior_samples(b4.1)
cov(post[, 1:2])

             b_Intercept        sigma
b_Intercept  0.164413984 -0.002461927
sigma       -0.002461927  0.090340147

Lets extract the HDPI

library(tidybayes)

post %>% 
  select(-lp__) %>% 
  gather(parameter) %>% 
  group_by(parameter) %>% 
  median_hdi(value)

Now lets add a predictor.

d3 <- d2 %>% 
  mutate(weight = (weight - mean(weight))/sd(weight)) %>% 
  select(height, weight)

b4.3 <- 
  brm(data = d3, family = gaussian,
      height ~ 1 + weight,
      prior = c(prior(normal(156, 100), class = Intercept),
                prior(normal(0, 10), class = b),
                prior(uniform(0, 50), class = sigma)),
      iter = 46000, warmup = 45000, chains = 4, cores = 4,
      control = list(adapt_delta = 0.8, 
                     max_treedepth = 10))

starting worker pid=46313 on localhost:11970 at 15:00:14.223
starting worker pid=46323 on localhost:11970 at 15:00:14.471
starting worker pid=46333 on localhost:11970 at 15:00:14.704
starting worker pid=46343 on localhost:11970 at 15:00:14.934

SAMPLING FOR MODEL 'a71d137dd5c6103d8c66ecffb73438c3' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 5.6e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.56 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:     1 / 46000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'a71d137dd5c6103d8c66ecffb73438c3' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 5.7e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.57 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:     1 / 46000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'a71d137dd5c6103d8c66ecffb73438c3' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 7.3e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.73 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:     1 / 46000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'a71d137dd5c6103d8c66ecffb73438c3' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 5e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.5 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:     1 / 46000 [  0%]  (Warmup)
Chain 1: Iteration:  4600 / 46000 [ 10%]  (Warmup)
Chain 2: Iteration:  4600 / 46000 [ 10%]  (Warmup)
Chain 3: Iteration:  4600 / 46000 [ 10%]  (Warmup)
Chain 4: Iteration:  4600 / 46000 [ 10%]  (Warmup)
Chain 1: Iteration:  9200 / 46000 [ 20%]  (Warmup)
Chain 2: Iteration:  9200 / 46000 [ 20%]  (Warmup)
Chain 3: Iteration:  9200 / 46000 [ 20%]  (Warmup)
Chain 1: Iteration: 13800 / 46000 [ 30%]  (Warmup)
Chain 4: Iteration:  9200 / 46000 [ 20%]  (Warmup)
Chain 2: Iteration: 13800 / 46000 [ 30%]  (Warmup)
Chain 3: Iteration: 13800 / 46000 [ 30%]  (Warmup)
Chain 4: Iteration: 13800 / 46000 [ 30%]  (Warmup)
Chain 1: Iteration: 18400 / 46000 [ 40%]  (Warmup)
Chain 2: Iteration: 18400 / 46000 [ 40%]  (Warmup)
Chain 3: Iteration: 18400 / 46000 [ 40%]  (Warmup)
Chain 3: Iteration: 23000 / 46000 [ 50%]  (Warmup)
Chain 2: Iteration: 23000 / 46000 [ 50%]  (Warmup)
Chain 1: Iteration: 23000 / 46000 [ 50%]  (Warmup)
Chain 4: Iteration: 18400 / 46000 [ 40%]  (Warmup)
Chain 3: Iteration: 27600 / 46000 [ 60%]  (Warmup)
Chain 2: Iteration: 27600 / 46000 [ 60%]  (Warmup)
Chain 1: Iteration: 27600 / 46000 [ 60%]  (Warmup)
Chain 3: Iteration: 32200 / 46000 [ 70%]  (Warmup)
Chain 4: Iteration: 23000 / 46000 [ 50%]  (Warmup)
Chain 2: Iteration: 32200 / 46000 [ 70%]  (Warmup)
Chain 1: Iteration: 32200 / 46000 [ 70%]  (Warmup)
Chain 3: Iteration: 36800 / 46000 [ 80%]  (Warmup)
Chain 2: Iteration: 36800 / 46000 [ 80%]  (Warmup)
Chain 1: Iteration: 36800 / 46000 [ 80%]  (Warmup)
Chain 4: Iteration: 27600 / 46000 [ 60%]  (Warmup)
Chain 3: Iteration: 41400 / 46000 [ 90%]  (Warmup)
Chain 2: Iteration: 41400 / 46000 [ 90%]  (Warmup)
Chain 1: Iteration: 41400 / 46000 [ 90%]  (Warmup)
Chain 3: Iteration: 45001 / 46000 [ 97%]  (Sampling)
Chain 2: Iteration: 45001 / 46000 [ 97%]  (Sampling)
Chain 3: Iteration: 46000 / 46000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 15.0375 seconds (Warm-up)
Chain 3:                0.54567 seconds (Sampling)
Chain 3:                15.5832 seconds (Total)
Chain 3: 
Chain 1: Iteration: 45001 / 46000 [ 97%]  (Sampling)
Chain 2: Iteration: 46000 / 46000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 15.7736 seconds (Warm-up)
Chain 2:                0.483538 seconds (Sampling)
Chain 2:                16.2572 seconds (Total)
Chain 2: 
Chain 4: Iteration: 32200 / 46000 [ 70%]  (Warmup)
Chain 1: Iteration: 46000 / 46000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 16.4164 seconds (Warm-up)
Chain 1:                0.604877 seconds (Sampling)
Chain 1:                17.0213 seconds (Total)
Chain 1: 
Chain 4: Iteration: 36800 / 46000 [ 80%]  (Warmup)
Chain 4: Iteration: 41400 / 46000 [ 90%]  (Warmup)
Chain 4: Iteration: 45001 / 46000 [ 97%]  (Sampling)
Chain 4: Iteration: 46000 / 46000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 24.0847 seconds (Warm-up)
Chain 4:                0.489564 seconds (Sampling)
Chain 4:                24.5742 seconds (Total)
Chain 4:

plot(b4.3)

Examine covariance

pairs(b4.3)

Creating predictions

post <- posterior_samples(b4.3)

mu_at_50 <-
  post %>% 
  transmute(mu_at_50 = b_Intercept + b_weight * 50)

mu_at_50 %>%
  ggplot(aes(x = mu_at_50)) +
  geom_density(size = 0, fill = "grey75") +
  stat_pointintervalh(aes(y = 0), 
                      point_interval = mode_hdi, .width = .95) +
  scale_y_continuous(NULL, breaks = NULL) +
  labs(x = expression(mu["height | weight = 50"])) +
  theme_classic()

Now we extract complete predictions including the mean and added variance.

weight_seq <- tibble(weight = seq(from = 25, to = 70, by = 1))

pred_height <-
  predict(b4.3,
          newdata = weight_seq) %>%
  as_tibble() %>%
  bind_cols(weight_seq)
  
pred_height %>%
  slice(1:6)

ggplot(pred_height, aes(weight, Estimate)) +
  geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), fill = "grey75") +
  geom_line() +
  theme(panel.grid = element_blank())

We can overlay this with the mean interval too

mu_summary <-
  fitted(b4.3, 
         newdata = weight_seq) %>%
  as_tibble() %>%
  bind_cols(weight_seq)

ggplot(mu_summary, aes(weight, Estimate)) +
  geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), fill = "grey75") +
  geom_line() +
  theme(panel.grid = element_blank())

inner_join(pred_height, mu_summary, by = "weight", suffix=c("_height", "_mean")) %>% str

Classes ‘tbl_df’, ‘tbl’ and 'data.frame':   46 obs. of  9 variables:
 $ Estimate_height : num  300 306 312 318 324 ...
 $ Est.Error_height: num  8.38 8.73 8.98 9.06 9.32 ...
 $ Q2.5_height     : num  285 289 295 300 305 ...
 $ Q97.5_height    : num  317 323 329 336 342 ...
 $ weight          : num  25 26 27 28 29 30 31 32 33 34 ...
 $ Estimate_mean   : num  300 306 312 318 324 ...
 $ Est.Error_mean  : num  6.68 6.95 7.21 7.48 7.75 ...
 $ Q2.5_mean       : num  288 293 298 304 309 ...
 $ Q97.5_mean      : num  313 320 326 332 339 ...

inner_join(pred_height, mu_summary, by = "weight", suffix=c("_height", "_mean")) %>% 
  ggplot(aes(x = weight, y = Estimate_height)) +
  geom_ribbon(aes(ymin = Q2.5_height, ymax = Q97.5_height), fill = "grey25") +
  geom_ribbon(aes(ymin = Q2.5_mean, ymax = Q97.5_mean), fill = "grey45") +
  geom_line() +
  theme(panel.grid = element_blank())

LS0tCnRpdGxlOiAiQ2hhcHRlciA0IC0gTGluZWFyIE1vZGVscyIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKYGBge3IgbWVzc2FnZT1GQUxTRSwgd2FybmluZz1GQUxTRX0KbGlicmFyeSh0aWR5dmVyc2UpCmxpYnJhcnkoc2tpbXIpCmxpYnJhcnkocmV0aGlua2luZykKYGBgCgojIE5vdGVzCgojIyBJbnRyb2R1Y3Rpb24KTGluZWFyIHJlZ3Jlc3Npb24gZGVzY3JpYmVzIGEgbWVhc3VyZW1lbnRzIG1lYW4gYW5kIHZhcmlhbmNlIGFzIHRoZSBhZGRpdGlvbiBvZiBvdGhlciBtZWFzdXJlbWVudHMuIEl0IGFzc3VtZXMgdGhlIGVycm9ycyBpbiB0aGUgcHJpbWFyeSBtZWFzdXJlbWVudCBhcmUgb2YgdGhlIGdhdXNzaWFuIGRpc3RyaWJ1dGlvbi4KCiMjIE5vcm1hbGl0eQoKMTAwMCBpbmRpdmlkdWFscyBjb25kdWN0IGEgcmFuZG9tIHdhbGsgb24gZWl0aGVyIHNpZGUgb2YgdGhlIGZvb3RiYWxsIGZpZWxkIHN0YXJ0aW5nIGF0IHRoZSA1MCB5YXJkIGxpbmUgKHZhbHVlID0gMCkuIFRoaXMgc2lnbmlmaWVzIHRoZSBub3JtYWwgZGlzdHJpYnV0aW9uLgpgYGB7ciBmaWcud2lkdGg9MTR9CklORElWSURVQUxTIDwtIDEwMDAKVFJJQUxTIDwtIDEwMAoKcHVycnI6Om1hcCgxOklORElWSURVQUxTLCB+IHJ1bmlmKFRSSUFMUywgLTEsIDEpKSAlPiUgCiAgcHVycnI6Om1hcChjdW1zdW0pICU+JSAKICByZWR1Y2UocmJpbmQpICU+JSAKICBhc190aWJibGUgJT4lIAogIGBjb2xuYW1lczwtYChwYXN0ZTAoMTpUUklBTFMsICJ0IikpICU+JSAKICBtdXRhdGUoaW5kaXZpZHVhbCA9IHJvd19udW1iZXIoKSkgJT4lIAogIHJlc2hhcGUyOjptZWx0KCJpbmRpdmlkdWFsIikgJT4lIAogIGFycmFuZ2UoaW5kaXZpZHVhbCwgdmFyaWFibGUpICU+JSAKICBnZ3Bsb3QoYWVzKHg9dmFyaWFibGUsIHk9dmFsdWUsIGdyb3VwPWluZGl2aWR1YWwpKSArIAogICAgZ2VvbV9saW5lKGFscGhhID0gSSgxLzIwKSwgc2l6ZT0xKSArCiAgICBnZW9tX2JveHBsb3QoYWVzKHg9dmFyaWFibGUsIHkgPSB2YWx1ZSksIGluaGVyaXQuYWVzID0gRkFMU0UsIGNvbG91ciA9ICJ5ZWxsb3ciLCBhbHBoYT0gSSgxLzQwKSkKYGBgCgpUaGUgYWJvdmUgZXhhbXBsZSBzaG93cyB1cyB0aGF0IG5vcm1hbGl0eSBpcyBhZGRpdGl2ZS4gCmBgYHtyfQojIEdlbmVyYXRlIDEyIG51bWJlcnMgYmV0d2VlbiAxIGFuZCAxLjEgdGhlbiBhZGQgdGhlbSBpbiB0byBvbmUgc3VtCiMgRG8gdGhpcyAxMDAwIHRpbWVzIGFuZCBjb2xsZWN0IHRoZSBudW1iZXJzIGluIGEgdmVjdG9yCmdyb3d0aCA8LSByZXBsaWNhdGUoMTAwMCwgc3VtKHJ1bmlmKDEyLDAsLjEpKSkKIyBQbG90IHRoZSBkaXN0cmlidXRpb24gb2YgdGhlcyBlbnVtYmVycy4gTm90ZSB0aGF0IGl0IGlzIG5vcm1hbC4KcXBsb3QoZ3Jvd3RoKQpgYGAKCldlIGNhbiBhbHNvIHNob3cgdGhhdCBub3JtYWxpdHkgaXMgbXVsdGlwbGljYXRpdmUuCmBgYHtyfQojIEdlbmVyYXRlIDEyIG51bWJlcnMgYmV0d2VlbiAxIGFuZCAxLjEgdGhlbiBtdWx0aXBseSB0aGVtIGluIHRvIG9uZSBwcm9kdWN0CiMgRG8gdGhpcyAxMDAwIHRpbWVzIGFuZCBjb2xsZWN0IHRoZSBudW1iZXJzIGluIGEgdmVjdG9yCmdyb3d0aCA8LSByZXBsaWNhdGUoMTAwMCwgcHJvZCgxICsgcnVuaWYoMTIsMCwuMSkpKQojcGxvdCB0aGUgZGlzdHJpYnV0aW9uIG9mIHRoZXNlIG51bWJlcnMuIE5vdGUgdGhhdCBpdCBpcyBub3JtYWwuCnFwbG90KElORElWSURVQUxTKQpgYGAKClRoaW5ncyBzdGFydCB0byBza2V3IGF0IHZlcnkgbGFyZ2UgZGV2aWF0ZXMsIHRob3VnaCB0aGlzIGNhbiBiZSBjb3JyZWN0ZWQgYnkgbG9nIHNjYWxpbmcuCmBgYHtyIH0KYmlnIDwtIHJlcGxpY2F0ZSgxMDAwMCwgcHJvZCgxICsgcnVuaWYoMTIsIDAsIC41KSkpCnNtYWxsIDwtIHJlcGxpY2F0ZSgxMDAwMCwgcHJvZCgxICsgcnVuaWYoMTIsIDAsIDAuMDEpKSkKCmRhdGFfZnJhbWUoYmlnID0gYmlnLCAKICAgICAgICAgICBzbWFsbCA9IHNtYWxsLCAKICAgICAgICAgICBsb2dfYmlnID0gbG9nKGJpZykpICU+JSAKICBnYXRoZXIgJT4lIAogIGdncGxvdChhZXModmFsdWUpKSArIAogICAgZ2VvbV9oaXN0b2dyYW0oKSArIAogICAgZmFjZXRfd3JhcCgia2V5Iiwgc2NhbGVzID0gImZyZWUiKQpgYGAKCiMjIEEgbGFuZ3VhZ2UgZm9yIGRlc2NyaWJpbmcgbW9kZWxzCgpFbGVtZW50czoKCjEuICoqQW4gb3V0Y29tZSB2YXJpYWJsZSoqOiBUaGUgc2V0IG9mIG1lYXN1cmVtZW50cyB3ZSBob3BlIHRvIHByZWRpY3Qgb3IgdW5kZXJzdGFuZAoyLiAqKkxpa2VsaWVob29kIGRpc3RyaWJ1dGlvbioqOiBGb3IgZWFjaCBvZiB0aGVzZSBvdXRjb21lIHZhcmlhYmxlcywgd2UgZGVmaW5lIHRoZSBwbGF1c2liaWxpdHkgb2YgdGhhdCBvYnNlcnZhdGlvbnMuCjMuICoqUHJlZGljdG9yIHZhcmlhYmxlcyoqOiBTZXQgb2Ygb3RoZXIgbWVhc3VyZW1lbnRzIHRoYXQgd2UgaG9wZSB0byB1c2UgdG8gcHJlZGljdCBvciB1bmRlcnN0YW5kIHRoZSBvdXRjb21lLgoKQXBwcm9hY2g6CgoxLiBSZWxhdGUgdGhlIGV4YWN0IHNoYXBlIG9mIHRoZSBsaWtlbGlob29kIGRpc3RyaWJ1dGlvbiB0byB0aGUgcHJlZGljdG9yIHZhcmlhYmxlcyB2aWEgcGFyYW1ldGVycwoyLiBDaG9vc2UgcHJpb3JzIGZvciBhbGwgdGhlIHBhcmFtZXRlcnMgaW4gdGhlIG1vZGVsLgoKCiMjIEV4YW1wbGU6IEEgZ2F1c3NpYW4gbW9kZWwgb2YgaGVpZ2h0CgpXZSdsbCBidWlsZCBhIHJlZ3Jlc3Npb24gdG8gcHJlZGljdCBhIGh1bWFuJ3MgaGlnaHQsIHdoaWNoIHdlIHdpbGwgbW9kZWwgYXMgYSBHYXVzc2lhbiBkaXN0cmlidXRpb24gd2l0aCB0d28gcGFyYW1ldGVycywgbWVhbiBhbmQgc3RhbmRhcmQgZGV2aWF0aW9uLiBUaGVyZSBhcmUgYW4gaW5maW5pdGUgbnVtYmVyIG9mIHBvc3NpYmxlIEdhdXNzaWFuIGRpc3RyaWJ1dGlvbnMuIFdlIHdhbnQgb3VyIEJheWVzaWFuIG1hY2hpbmUgdG8gY29uc2lkZXIgZXZlcnkgcG9zc2libGUgZGlzdHJpYnV0aW9uLCBlYWNoIGRlZmluZWQgYnkgYSBjb21iaW5hdGlvbiBvZiBtZWFuZCBhbmQgc2QsIGFuZCByYW5rIHRoZW0gYnkgcG9zdGVyaW9yIHBsYXVzaWJpbGl0eSBnaXZlbiB0aGUgZGF0YS4KCmBgYHtyfQpkYXRhKEhvd2VsbDEpCmQgPC0gSG93ZWxsMQpkMiA8LSBkICU+JSAKICBmaWx0ZXIoYWdlID49IDE4KQpkMiAlPiUgCiAgc2tpbQpgYGAKCkxldHMgc3BlY2lmaXkgb3VyIG1vZGVsIGxpa2Ugc28uIFRoZSBmaXJzdCBsaW5lIGlzIHRoZSBsaWtlbGlob29kIG9mIG91ciBvdXRjb21lLCBoZWlnaHQuIFRoZSBzZWNvbmQgdHdvIGxpbmVzIGFyZSBwcmlvcnMgZm9yIHRoZSBwYXJhbWV0ZXJzIGluIG91ciBsaWtlbGlob29kLgokJApoX2kgXHNpbSBOb3JtYWwoXG11LCBcc2lnbWEpIFxcClxtdSBcc2ltIE5vcm1hbCgxNzgsIDIwKSBcXApcc2lnbWEgXHNpbSBVbmlmb3JtKDAsIDUwKQokJAoKQXMgaGVpZ2h0IGlzIGEgcGh5c2ljYWwgcGhlbm9tZW5vbiwgd2UgY2FuIHVzZSBvdXIgaW50dWl0aW9uIHRvIGFzc2lnbiBhIHByaW9yLgpgYGB7cn0KZGF0YV9mcmFtZSh4PTEwMDoyNTAsIAogICAgICAgICAgIHk9cHVycnI6Om1hcF9kYmwoMTAwOjI1MCwgfiBkbm9ybSgueCwgMTc4LCAyMCkpKSAlPiUgCiAgZ2dwbG90KGFlcyh4LHkpKSArIGdlb21fcGF0aCgpCmBgYAoKV2UgZ2l2ZSB0aGUgc3RhbmRhcmQgZGV2aWF0aW9uIGEgdHJ1bHkgZmxhdCBwcmlvci4gV2UgcmVhbGx5IGp1c3Qgd2FudCB0byBtYWtlIHN1cmUgaXRzIHBvc2l0aXZlLgpgYGB7cn0KZGF0YV9mcmFtZSh4PS0xMDo2MCwgCiAgICAgICAgICAgeT1wdXJycjo6bWFwX2RibCgtMTA6NjAsIH4gZHVuaWYoLngsIDAsIDUwKSkpICU+JSAKICBnZ3Bsb3QoYWVzKHgseSkpICsgZ2VvbV9wYXRoKCkKYGBgCgpTcGVjaWZ5aW5nIGEgcHJpb3IgZm9yIHRoZSBwYXJhbWV0ZXJzIG9mIG91ciBvdXRjb21lIGltcGxpY2l0bHkgZ2l2ZSB1cyBhIHByaW9yIGZvciBvdXIgb3V0Y29tZS4KYGBge3J9CnNhbXBsZV9tdSA8LSBybm9ybSgxZTQsIDE3OCwgMjApCnNhbXBsZV9zaWdtYSA8LSBydW5pZigxZTQsIDAsIDUwKQpwcmlvcl9oIDwtIHJub3JtKDFlNCwgc2FtcGxlX211LCBzYW1wbGVfc2lnbWEpCnFwbG90KHByaW9yX2gpCmBgYAoKTGV0cyB1cyBncmlkIGFwcHJveGltYXRpb24gdG8gZ2VuZXJhdGUgdGhlIHBvc3Rlcmlvci4gVG8gZG8gc28sIHdlIGVudW1lcmF0ZSBhIGxpc3Qgb2YgY2FuZGlkYXRlIHZhbHVlcyBmb3Igb3VyIHBhcmFtZXRlcnMgYW5kIHN0b3JlIHRoZW0gaW4gcG9zdC4gV2UgdGhlbiBjYWxjdWxhdGUgdGhlIGxvZyBsaWtlbGlob29kIGZvciBlYWNoIG9mIHRob3NlIHBhcmFtZXRlcnMgZ2l2ZW4gYWxsIHRoZSBoZWlnaHRzIHdlIGhhdmUgc2VlbiBhbmQgc3RvcmUgdGhlbSBpbiB0aGUgY29sdW1uIExMLiBXZSB1c2UgbG9nIHRyYW5zZm9ybWF0aW9uIHRvIGF2b2lkIHJvdW5kaW5nIGVycm9ycy4gVGhlbiB3ZSBtdWx0aXBseSBieSB0aGUgcHJpb3JzIChsb2cgYWRkaXRpb24gaXMgdGhlIHNhbWUgYXMgbXVsdGlwbHlpbmcpCmBgYHtyfQojIEdlbmVyYXRlIGNhbmRpZGF0ZSBwYXJhbWV0ZXJzCm11X2xpc3QgPC0gc2VxKGZyb209MTQwLCB0bz0xNjAsIGxlbmd0aC5vdXQ9MjAwKQpzaWdtYV9saXN0IDwtIHNlcShmcm9tPTQsIHRvPTksIGxlbmd0aC5vdXQ9MjAwKQpwb3N0IDwtIGV4cGFuZC5ncmlkKG11PW11X2xpc3QsIHNpZ21hPXNpZ21hX2xpc3QpCgojIENhbGN1bGF0ZSB0aGUgbGlrZWxpaG9vZCBvZiBldmVyeSBwb3NzaWJsZQojIGNvbWJpbmF0aW9uIGNhbmRpZGF0ZSBwYXJhbWV0ZXJzCnBvc3QkTEwgPC0gc2FwcGx5KDE6bnJvdyhwb3N0KSwgZnVuY3Rpb24oaSkgewogIHN1bShkbm9ybSgKICAgIGQyJGhlaWdodCwKICAgIG1lYW49cG9zdCRtdVtpXSwKICAgIHNkPXBvc3Qkc2lnbWFbaV0sCiAgICBsb2c9VFJVRSkpCiAgfSkKCiMgTXVsdGlwbHkgcHJpb3JzCnBvc3QkcHJvZCA8LSBwb3N0JExMICsgZG5vcm0ocG9zdCRtdSwgMTc4LCAyMCwgVFJVRSkgKwogIGR1bmlmKHBvc3Qkc2lnbWEsIDAsIDUwLCBUUlVFKQoKIyBDYWxjdWxhdGVkIHNjYWxlZCBwb3N0ZXJpb3IKcG9zdCRwcm9iIDwtIGV4cChwb3N0JHByb2QgLSBtYXgocG9zdCRwcm9kKSkKCiMgQ29udG91ciBwbG90IG9mIHByb2JhYmlsaXRpZXMKZ2dwbG90KHBvc3QsIGFlcyh4PW11LCB5PXNpZ21hLCB6PXByb2IpKSArIAogIGdlb21fY29udG91cihhZXMoY29sb3VyPXN0YXQobGV2ZWwpKSkgKyAKICBjb29yZF9jYXJ0ZXNpYW4oeGxpbT1jKDE0MCwgMTYwKSwgeWxpbT1jKDQsIDkpKQpgYGAKCk5vdyBsZXRzIHNhbXBsZSBmcm9tIHRoZSBwb3N0ZXJpb3IuCmBgYHtyfQojIFNhbXBsZSBjYW5kaWRhdGUgdmFsdWVzIGJ5IHBvc3RlcmlvcgpzYW1wbGVfcm93cyA8LSBzYW1wbGUoMTpucm93KHBvc3QpLCBzaXplPTFlNCwgcmVwbGFjZT1UUlVFLCBwcm9iPXBvc3QkcHJvYikKCiMgUGxvdCB0aGVpciAyZCBoaXN0b2dyYW0KZGF0YV9mcmFtZSgKICBzYW1wbGVfbXUgPSBwb3N0JG11W3NhbXBsZV9yb3dzXSwKICBzYW1wbGVfc2lnbWEgPSBwb3N0JHNpZ21hW3NhbXBsZV9yb3dzXQopICU+JSAKICBnZ3Bsb3QoYWVzKHNhbXBsZV9tdSwgc2FtcGxlX3NpZ21hKSkgKwogICAgZ2VvbV9iaW4yZChiaW5zPTMwKSArIAogICAgY29vcmRfY2FydGVzaWFuKHhsaW09YygxNDAsIDE2MCksIHlsaW09Yyg0LCA5KSkKYGBgCgpMZXRzIGFsc28gZXhhbWluZSB0aGUgbWFyZ2lucy4KYGBge3J9CnJlcXVpcmUoZ3JpZEV4dHJhKQoKcDEgPC0gcG9zdCRtdVtzYW1wbGVfcm93c10gJT4lIAogIHFwbG90ICsgZ2d0aXRsZSgibXUiKQpwMiA8LSBwb3N0JHNpZ21hW3NhbXBsZV9yb3dzXSAlPiUgCiAgcXBsb3QgKyBnZ3RpdGxlKCJzaWdtYSIpCgpncmlkLmFycmFuZ2UocDEsIHAyLCBuY29sPTIpCmBgYAoKYG1hcGAgZnJvbSB0aGUgcmV0aGlua2luZyBwYWNrYWdlIGFsbG93cyB1cyB0byBzb2x2ZSBmb3IgdGhlIHBvc3RlcmlvciB1c2luZyBxdWFkcmF0aWMgYXBwcm94aW1hdGlvbiBpbnN0ZWFkIG9mIGdyaWQgYXBwcm94aW1hdGlvbnMuIEl0IHNob3dzIHVzIGdhdXNzaWFuIGFwcHJveGltYXRpb25zIGZvciBlYWNoIHBhcmFtZXRlcidzIG1hcmdpbmFsIGRpc3RyaWJ1dGlvbi4KYGBge3J9CmZsaXN0IDwtIGFsaXN0KAogIGhlaWdodCB+IGRub3JtKG11LCBzaWdtYSksCiAgbXUgfiBkbm9ybSgxNzgsIDIwKSwKICBzaWdtYSB+IGR1bmlmKDAsIDUwKQopCm00LjEgPC0gcmV0aGlua2luZzo6bWFwKGZsaXN0LCBkYXRhPWQyKQpwcmVjaXMobTQuMSkKYGBgCgpWYXJpYW5jZSBjb3ZhcmlhbmNlIG1hdHJpeCB0ZWxscyB1cyBob3cgZWFjaCBwYXJhbWV0ZXIgcmVsYXRlcyB0byBldmVyeSBvdGhlciBwYXJhbWV0ZXIgaW4gdGhlIHBvc3RlcmlvciBkaXN0cmlidXRpb24uCmBgYHtyfQp2Y292KG00LjEpCiMgVGhpcyBjYW4gYmUgZmFjdG9yZWQgaW4gdG8gYSB2ZWN0b3Igb2YgdmFyaWFuY2VzCmRpYWcodmNvdihtNC4xKSkKIyBPciBhIGNvcnJlbGF0aW9uIG1hdHJpeApjb3YyY29yKHZjb3YobTQuMSkpCmBgYAoKTGV0cyBkcmF3IGZyb20gdGhlIHBvc3RlcmlvciBieSBzYW1wbGluZyB2ZWN0b3JzIG9mIHZhbHVlcyBmcm9tIGEgbXVsdGktZGltZW5zaW9uYWwgZ2F1c3NpYW4gZGlzdHJpYnV0aW9uLgpgYGB7cn0KcG9zdCA8LSBleHRyYWN0LnNhbXBsZXMobTQuMSwgbj0xZTQpCnByZWNpcyhwb3N0KQpgYGAKCiMjIyBBZGRpbmcgYSBwcmVkaWN0b3IgdmFyaWFibGUsIHdlaWdodAoKSXQgbG9va3MgbGlrZSB3ZWlnaHQgd2lsbCBiZSBhIGdvb2QgcHJlZGljdG9yIGZvciBoZWlnaHQgaW4gYWR1bHRzLgpgYGB7cn0KZ2dwbG90KGQyLCBhZXMoeD13ZWlnaHQsIHk9aGVpZ2h0KSkgKyAKICBnZW9tX3BvaW50KHNoYXBlPTEpICsKICBnZW9tX3Ntb290aChtZXRob2Q9bG0pCmBgYAoKVGhlIGdvYWwgaGVyZSBpcyB0byBtYWtlIHRoZSBwYXJhbWV0ZXIgZm9yIHRoZSBtZWFuIG9mIGEgR2F1c3NpYW4gZGlzdHJpYnV0aW9uLCAkXG11JCwgaW50byBhIGxpbmVhciBmdW5jdGlvbiBvZiB0aGUgcHJlZGljdG9yIHZhcmlhYmxlLiBUaGlzIGVuY29kZXMgdGhlIGFzc3VtcHRpb24gdGhhdCB0aGUgcHJlZGljdG9yIHZhcmlhYmxlIGhhcyBhIHBlcmZlY3RseSBjb25zdGFudCBhbmQgYWRkaXRpdmUgcmVsYXRpb25zaGlwIHRvIHRoZSBtZWFuIG9mIHRoZSBvdXRjb21lLgpgYGB7cn0KbTQuMyA8LSByZXRoaW5raW5nOjptYXAoCiAgYWxpc3QoCiAgICBoZWlnaHQgfiBkbm9ybShtdSwgc2lnbWEpLAogICAgbXUgPC0gYSArIGIqd2VpZ2h0LAogICAgYSB+IGRub3JtKDE1NiwgMTAwKSwKICAgIGIgfiBkbm9ybSgwLCAxMCksCiAgICBzaWdtYSB+IGR1bmlmKDAsIDUwKSksCiAgZGF0YSA9IGQyKQojIERpc3BsYXkgdGFibGUgb2YgZXN0aW1hdGVzCnByZWNpcyhtNC4zLCBjb3JyPVRSVUUpCmBgYAoKVGhlIHRhYmxlIGFib3ZlIHRlbGxzIHVzIHRoYXQgYSBwZXJzb24gMSBrZyBoZWF2aWVyIGlzIGV4cGVjdGVkIHRvIGJlIC45IGNtIHRhbGxlci4gSXQgYWxzbyB0ZWxscyB1cyB0aGF0IGEgcGVyc29uIHdpdGggMCB3ZWlnaHQgaXMgMTE0IGNtIHRhbGwsIHdoaWNoIGlzIG9mIGNvdXJzZSBub25zZW5zZTsgVGhpcyBpcyB3aHkgaXQgaXMgdXN1YWxseSBpbXBvcnRhbnQgdG8gaGF2ZSB2ZXJ5IHdlYWsgcHJpb3JzIGZvciBpbnRlcmNlcHRzIGluIG1hbnkgY2FzZXMuIEZpbmFsbHksIGl0IGV4cGxhaW5zIHRoYXQgOTUlIG9mIHBsYXVzaWJsZSBoZWlnaHRzIGxpZSB3aXRoaW4gMTAgY20gb2YgdGhlIG1lYW4uCgpDZW50ZXJpbmcgaXMgdGhlIHByb2NlZHVyZSBvZiBzdWJ0cmFjdGluZyB0aGUgbWVhbiBvZiBhIHZhcmlhYmxlIGZyb20gZWFjaCB2YWx1ZS4gCmBgYHtyfQpkMiR3ZWlnaHQuYyA8LSBkMiR3ZWlnaHQgLSBtZWFuKGQyJHdlaWdodCkKCm00LjQgPC0gcmV0aGlua2luZzo6bWFwKAogIGFsaXN0KAogICAgaGVpZ2h0IH4gZG5vcm0obXUsIHNpZ21hKSwKICAgIG11IDwtIGEgKyBiKndlaWdodC5jLAogICAgYSB+IGRub3JtKDE3OCwgMTAwKSwKICAgIGIgfiBkbm9ybSgwLCAxMCksCiAgICBzaWdtYSB+IGR1bmlmKDAsIDUwKSksCiAgZGF0YSA9IGQyKQojIERpc3BsYXkgdGFibGUgb2YgZXN0aW1hdGVzCnByZWNpcyhtNC40LCBjb3JyPVRSVUUpCmBgYAoKQnkgKipjZW50ZXJpbmcqKiBvdXIgcHJlZGljdG9yIHZhcmlhYmxlLCB0aGUgaW50ZXJwcmV0YXRpb24gb2YgdGhlIGludGVyY2VwdCBoYXMgbm93IGJlY29tZSB0aGUgZXhwZWN0ZWQgdmFsdWUgb2YgdGhlIG91dGNvbWUgd2hlbiB0aGUgcHJlZGljdG9yIGlzIGF0IGl0cyBhdmVyYWdlIHZhbHVlLiBUaGlzIG1ha2VzIGludGVycHJldGluZyB0aGUgaW50ZXJjZXB0IGEgbG90IGVhc2llci4gCgojIyMgVmlzdWFsaXppbmcgb3VyIGVzdGltYXRlcwoKV2UgY2FuIHBsb3Qgb3VyIG1heGltdW0gYSBwb3N0ZXJpb3JpIGZpdCBoZXJlCgpgYGB7cn0KcCA8LSBnZ3Bsb3QoZDIsIGFlcyh4PXdlaWdodCwgeT1oZWlnaHQpKSArCiAgZ2VvbV9wb2ludChzaGFwZT0xLCBzaXplPTIpICsKICBnZW9tX2FibGluZShhZXMoaW50ZXJjZXB0ID0gY29lZihtNC4zKVsiYSJdLAogICAgICAgICAgICAgICAgICBzbG9wZSA9IGNvZWYobTQuMylbImIiXSkpCnAKYGBgCgpUaGlzIGxpbmUgb25seSByZXByZXNlbnRzIG9uZSBwaWVjZSBvZiB0aGUgcG9zdGVyaW9yLiBXZSBjb3VsZCBpbnN0ZWFkIHNhbXBsZSB2YWx1ZXMgb2YgJFxhbHBoYSQgYW5kICRcYmV0YSQgZnJvbSB0aGUgcG9zdGVyaW9yIHRvIGdlbmVyYXRlIG1hbnkgbGluZXMgdG8gYnVpbGQgYm91bmRzIG9mIHVuY2VydGFpbnR5LiBXZSdsbCBkbyB0aGlzIHdpdGggYSBtdWNoIHNtYWxsZXIgZGF0YXNldCB0byBtYWduaWZ5IGVmZmVjdHMuCmBgYHtyfQpOIDwtIDEwCmROIDwtIGQyWyAxOk4gLCBdCm1OIDwtIHJldGhpbmtpbmc6Om1hcCgKICBhbGlzdCgKICAgIGhlaWdodCB+IGRub3JtKCBtdSAsIHNpZ21hICkgLAogICAgbXUgPC0gYSArIGIqd2VpZ2h0ICwKICAgIGEgfiBkbm9ybSggMTc4ICwgMTAwICkgLAogICAgYiB+IGRub3JtKCAwICwgMTAgKSAsCiAgICBzaWdtYSB+IGR1bmlmKCAwICwgNTAgKQopICwgZGF0YT1kTiApCgpwb3N0IDwtIGV4dHJhY3Quc2FtcGxlcyhtTiwgbj0yMCkKaGVhZChwb3N0KQpgYGAKCmBgYHtyfQpwIDwtIGdncGxvdChkTiwgYWVzKHdlaWdodCwgaGVpZ2h0KSkgKwogIGdlb21fcG9pbnQoc2hhcGU9MSkKCmZvciAoaSBpbiAxOjIwKSB7CiAgcCA8LSBwICsgZ2VvbV9hYmxpbmUoaW50ZXJjZXB0ID0gcG9zdCRhW2ldLCBzbG9wZSA9IHBvc3QkYltpXSwgYWxwaGEgPSBJKDMvMTApKQp9CgpwCmBgYAoKVG8gZG8gdGhpcyBtb3JlIGdlbmVyYWxseToKCjEuIFNhbXBsZSBwYXJhbWV0ZXIgdmFsdWVzIGZyb20gdGhlIHBvc3Rlcmlvci4KMi4gVXNlIHRoZSBwYXJhbWV0ZXIgdmFsdWVzIHRvIGNyZWF0ZSBwcmVkaWN0aW9ucyBvZiB0aGUgcGFyYW1ldGVycyBvZiBpbnRlcmVzdCwgaW5jbHVkaW5nIHRoZSBmaW5hbCBvdXRjb21lLgozLiBVc2Ugc3VtbWFyeSBmdW5jdGlvbnMgbGlrZSBtZWFuLCBIRE9QLCBQSSB0byBmaW5kIGF2ZXJhZ2VzIGFuZCBsb3dlciBhbmQgdXBwZXIgYm91bmRzIG9mIG91ciB1bmNlcnRhaW50eS4KCgpgYGB7cn0KIyBEcmF3IHBhcmFtZXRlciB2YWx1ZXMgZnJvbSB0aGUgcG9zdGVyaW9yIGRpc3RyaWJ1dGlvbgpwb3N0ZXJpb3Jfc2FtcGxlcyA8LSBleHRyYWN0LnNhbXBsZXMobTQuMykKCiMgQ3JlYXRlIG91dCB4LWF4aXMgb2YgaW50ZXJlc3QKd2VpZ2h0c19heGlzIDwtIDI1OjcwCgojIEZ1bmN0aW9uIHRvIGNhbGN1bGF0ZSBmaXQgdmFsdWVzIGZvciBcbXUKbXVfbGluayA8LSBmdW5jdGlvbih3ZWlnaHQpIHBvc3QkYSArIHBvc3QkYip3ZWlnaHQKbXUgPC0gcHVycnI6Om1hcCh3ZWlnaHRzX2F4aXMsIG11X2xpbmspCgojIEZ1bmN0aW9uIHRvIGNhbGN1bGF0ZSBmaXQgaGVpZ2h0cwojIFRoaXMgaW5jb3Jwb3JhdGVzIHZhcmlhYmlsaXR5IGFzIHdlbGwKc2ltdWxhdGUgPC0gZnVuY3Rpb24od2VpZ2h0LCBwYXJhbWV0ZXJfc2FtcGxlcykgewogIHJub3JtKG4gPSBucm93KHBhcmFtZXRlcl9zYW1wbGVzKSwKICAgICAgICBtZWFuID0gcGFyYW1ldGVyX3NhbXBsZXMkYSArIHBhcmFtZXRlcl9zYW1wbGVzJGIgKiB3ZWlnaHQsCiAgICAgICAgc2QgPSBwYXJhbWV0ZXJfc2FtcGxlcyRzaWdtYSl9CnNpbXVsYXRlZF9oZWlnaHRzIDwtIHB1cnJyOjptYXAod2VpZ2h0c19heGlzLCB+IHNpbXVsYXRlKC54LCBwb3N0ZXJpb3Jfc2FtcGxlcykpCgojIENvbWJpbmUgYWxsIHRoaXMgaW5mb3JtYXRpb24gdG8gcGxvdApkYXRhX2ZyYW1lKAogICMgWCBheGlzCiAgd2VpZ2h0ID0gd2VpZ2h0c19heGlzLAogICMgRXN0aW1hdGVzIGZvciBtZWFuIGhlaWdodAogIG11X21lYW4gPSBtYXBfZGJsKG11LCBtZWFuKSwKICBtdV9sb3dlciA9IG1hcF9kYmwobXUsIH5IUERJKC54LCAuODkpWzFdKSwKICBtdV91cHBlciA9IG1hcF9kYmwobXUsIH5IUERJKC54LCAuODkpWzJdKSwKICAjIEVzdGltYXRlcyBmb3IgaGVpZ2h0CiAgaGVpZ2h0X21lYW4gPSBtYXBfZGJsKHNpbXVsYXRlZF9oZWlnaHRzLCBtZWFuKSwKICBoZWlnaHRfbG93ID0gbWFwX2RibChzaW11bGF0ZWRfaGVpZ2h0cywgfkhQREkoLngsIC44OSlbMV0pLAogIGhlaWdodF9oaWdoID0gbWFwX2RibChzaW11bGF0ZWRfaGVpZ2h0cywgfkhQREkoLngsIC44OSlbMl0pKSAlPiUgCiAgZ2dwbG90KGRhdGE9LikgKwogIGdlb21fcG9pbnQoaW5oZXJpdC5hZXMgPSBGQUxTRSwgZGF0YSA9IGQyLAogICAgICAgICAgICAgYWVzKHg9d2VpZ2h0LCB5PWhlaWdodCksIHNoYXBlPTEpICsKICBnZW9tX2xpbmUoYWVzKHg9d2VpZ2h0LCB5PW11X21lYW4pKSArCiAgZ2VvbV9yaWJib24oYWVzKHggPSB3ZWlnaHQsIHltaW49bXVfbG93ZXIsIHltYXg9bXVfdXBwZXIpLCBhbHBoYT0uNSkgKwogIGdlb21fcmliYm9uKGFlcyh4ID0gd2VpZ2h0LCB5bWluID0gaGVpZ2h0X2xvdywgeW1heCA9IGhlaWdodF9oaWdoKSwgYWxwaGE9LjMpICsKICBnZ3RpdGxlKCJUaGUgZWZmZWN0IG9mIHdlaWdodCBvbiBoZWlnaHQiLCBzdWJ0aXRsZSA9ICJVc2luZyA4OSUgaGlnaGVzdCBwb3N0ZXJpb3IgZGVuc2l0eSBpbnRlcnZhbHMiKSArCiAgY29vcmRfY2FydGVzaWFuKHhsaW09YygzMCwgNjApLHlsaW09YygxMzUsIDE4MCkpCmBgYAoKIyMgUG9seW5vbWlhbCBSZWdyZXNzaW9uCgpTdGFuZGFyZGl6aW5nIG91ciBwcmVkaWN0b3JzIG1lYW5zIHNjYWxpbmcgYnV0IGFsc28gZGl2aWRpbmcgYnkgdGhlIHN0YW5kYXJkIGRldmlhdGlvbi4gVGhpcyBjaGFuZ2VzIG91ciBpbnRlcnByZXRhdGlvbiBvZiBjb250aW51b3VzIHZhcmlhYmxlcyB0byBtZWFuIHRoZSB1bml0IGNoYW5nZSBpbiBvdXRjb21lIGZyb20gYSBjYWhuZ2Ugb2Ygb25lIHN0YW5kYXJkIGRldmlhdGlvbi4gCgpgYGB7cn0KIyBTdGFuZGFyZGl6ZSBjb250aW51b3VzIHZhcmlhYmxlCnN0YW5kYXJkaXplIDwtIGZ1bmN0aW9uKHZlYykgewogIGNlbnRlcmVkIDwtIHZlYyAtIG1lYW4odmVjKQogIGNlbnRlcmVkIC8gc2QodmVjKQp9CmQkd2VpZ2h0LnMgPC0gc3RhbmRhcmRpemUoZCR3ZWlnaHQpCgojIERlZmluZSBhbm90aGVyIHBvbHlub21pYWwgdGVybQpkJHdlaWdodC5zMiA8LSBkJHdlaWdodC5zXjIKCm00LjUgPC0gbWFwKAogIGFsaXN0KAogICAgIyBPdXRjb21lCiAgICBoZWlnaHQgfiBkbm9ybShtdSwgc2lnbWEpLAogICAgIyBSZWdyZXNzaW9uIHNwZWNpZmljYXRpb24KICAgIG11IDwtIGEgKyBiMSp3ZWlnaHQucyArIGIyKndlaWdodC5zMiwKICAgIGEgfiBkbm9ybSgxNzggLCAxMDApLAogICAgIyBUaGVzZSBwcmlvcnMgYXJlIHZlcnkgd2VhayEKICAgIGIxIH4gZG5vcm0oMCAsIDEwKSwKICAgIGIyIH4gZG5vcm0oMCAsIDEwKSwKICAgIHNpZ21hIH4gZHVuaWYoMCAsIDUwKSksCiAgZGF0YT1kICkKcHJlY2lzKCBtNC41ICkKYGBgCgojIyBPbiB0cmFuc2Zvcm1hdGlvbnMKCkNlbnRlcmluZyBhbmQgc2NhbGluZyBhIGNvbnRpbnVvdXMgdmFyaWFibGUgaW5jcmVhc2VzIGludGVycHJldGFiaWxpdHkuIERlbWVhbmluZyBhIHByZWRpY3RvciB2YXJpYWJsZSBzdWJ0cmFjdHMgZXZlcnkgdmFsdWUgYnkgaXQncyBtZWFuLiBCeSBkZW1lYW5pbmcgYSBwcmVkaWN0b3IgdmFyaWFibGVzLCB0aGUgaW50ZXJjZXB0IHJlcHJlc2VudHMgdGhlIHZhbHVlIG9mIHRoZSBvdXRjb21lIHdoZW4gYWxsIHZhcmlhYmxlcyBhcmUgYXQgdGhlaXIgbWVhbi4gT25lIHN0YW5kYXJkaXplcyBhIHZhcmlhYmxlIGJ5IGRpdmlkaW5nIGV2ZXJ5IHZhbHVlIGJ5IHRoZSB2YXJpYWJsZXMgc3RhbmRhcmQgZGV2aWF0aW9uLiBTdGFuZGFyZGl6aW5nIGNhdXNlcyB0aGUgaW50ZXJwcmV0YXRpb24gb2YgY29lZmZpY2llbnRzIHRvIGJlIGluIHRlcm1zIG9mIHVuaXQgY2hhbmdlcyBpbiBTRCBvZiB0aGUgdmFyaWFibGUuCgoKTG9nIHRyYW5zZm9ybWluZyB0aGUgb3V0Y29tZSB2YXJpYWJsZSBtYWtlcyB5IG91ciBjb2VmZmljaWVudHMgcmVwcmVzZW50IHBlcmNlbnQgaW5jcmVhc2VzIGluIHRoZSBvdXRjb21lLiBMb2cgdHJhbnNmb3JtaW5nIGEgcHJlZGljdG9yIHZhcmlhYmxlIGRlc2NyaWJlcyBob3cgYSBzaW5nbGUgcGVyY2VudGFnZSBwb2ludCBpbmNyZWFzZSBpbiB0aGUgcHJlZGljdG9yIGFmZmVjdHMgeS4gTG9nIHRyYW5zZm9ybWluZyBib3RoIHRoZSBwcmVkaWN0b3IgYW5kIG91dGNvbWUgZGVzY3JpYmVzIGhvdyBwZXJjZW50YWdlIGNoYW5nZXMgaW4gb25lIGNyZWF0ZSBwZXJjZW50YWdlIGNoYW5nZXMgaW4gdGhlIG90aGVyLgoKCiMgSG9tZXdvcmsKCiMjIE1lZGl1bQojIyMgNE0xCkZvciB0aGUgbW9kZWwgZGVmaW5pdGlvbiBiZWxvdywgc2ltdWxhdGUgb2JzZXJ2ZWQgaGVpZ2h0cyBmcm9tIHRoZSBwcmlvci4KYGBge3J9Cm51bV9vYnNlcnZhdGlvbnMgPC0gMWU0Cgpwcmlvcl9tZWFuIDwtIHJub3JtKG51bV9vYnNlcnZhdGlvbnMsIDAsIDEwKQpwcmlvcl9zZCA8LSBydW5pZihudW1fb2JzZXJ2YXRpb25zLCAwLCAxMCkKCnByaW9yX2hlaWdodHMgPC0gcm5vcm0obnVtX29ic2VydmF0aW9ucywgcHJpb3JfbWVhbiwgcHJpb3Jfc2QpCgpxcGxvdChwcmlvcl9oZWlnaHRzKQpgYGAKCgojIyMgNE0yCmBgYHtyfQptb2RlbCA8LSBhbGlzdCgKICB5IH4gZG5vcm0obXUsIHNpZ21hKSwKICBtdSB+IGRub3JtKDAsIDEwKSwKICBzaWdtYSB+IGR1bmlmKDAsIDEwKQopCmBgYAoKCgojIyBIYXJkCiMjIyA0SDEKYGBge3J9CmRhdGEoIkhvd2VsbDEiKQpkIDwtIEhvd2VsbDEKZCR3ZWlnaHQgPC0gc3RhbmRhcmRpemUoZCR3ZWlnaHQpCmQKYGBgCgoKYGBge3J9Cm1vZGVsIDwtIG1hcCgKICBhbGlzdCgKICAgIGhlaWdodCB+IGRub3JtKG11LCBzaWdtYSksCiAgICBtdSA8LSBhICsgYiAqIHdlaWdodCwKICAgIGEgfiBkbm9ybShtZWFuID0gMCwgc2QgPSAxMDApLAogICAgYiB+IGRub3JtKG1lYW4gPSAwLCBzZCA9IDEwKSwKICAgIHNpZ21hIH4gZHVuaWYobWluID0gMCwgbWF4ID0gNjQpCiAgKSwKICBkYXRhID0gZAopCm1vZGVsCmBgYAoKYGBge3J9CiMgU3RhbmRhcmRpemUgb3VyIGlucHV0Cm5ld193ZWlnaHRzIDwtIHN0YW5kYXJkaXplKGMoNDYuOTUsIDQzLjcyLCA2NC43OCwgMzIuNTksIDU0LjYzKSkKCiMgRHJhdyBwYXJhbWV0ZXIgdmFsdWVzIGZyb20gdGhlIHBvc3RlcmlvciBkaXN0cmlidXRpb24KcG9zdGVyaW9yX3NhbXBsZXMgPC0gZXh0cmFjdC5zYW1wbGVzKG1vZGVsKQoKIyBGdW5jdGlvbiB0byBjYWxjdWxhdGUgZml0IGhlaWdodHMKc2ltdWxhdGVfaGVpZ2h0cyA8LSBmdW5jdGlvbih3ZWlnaHQsIHBhcmFtZXRlcl9zYW1wbGVzKSB7CiAgcm5vcm0oCiAgICBuID0gbnJvdyhwYXJhbWV0ZXJfc2FtcGxlcyksCiAgICBtZWFuID0gcGFyYW1ldGVyX3NhbXBsZXMkYSArIHBhcmFtZXRlcl9zYW1wbGVzJGIgKiB3ZWlnaHQsCiAgICBzZCA9IHBhcmFtZXRlcl9zYW1wbGVzJHNpZ21hKX0KCiMgR2VuZXJhdGUgYSBwcmVkaWN0aW9uIGZvciBlYWNoIG5ldyB3ZWlnaHQKIyBVc2UgZXZlcnkgcGFyYW1ldGVyIHNldCBpbiB0aGUgcG9zdGVyaW9yIHRvIGRvIHNvCnNpbXVsYXRlZF9oZWlnaHRzIDwtIHB1cnJyOjptYXAobmV3X3dlaWdodHMsIH4gc2ltdWxhdGUoLngsIHBvc3Rlcmlvcl9zYW1wbGVzKSkKCmRhdGFfZnJhbWUoCiAgaW5kaXZpZHVhbCA9IDE6NSwKICB3ZWlnaHQgPSBjKDQ2Ljk1LCA0My43MiwgNjQuNzgsIDMyLjU5LCA1NC42MyksCiAgZXhwZWN0ZWRfaGVpZ2h0ID0gbWFwX2RibChzaW11bGF0ZWRfaGVpZ2h0cywgbWVhbiksCiAgbG93ZXJfcGkgPSBtYXBfZGJsKHNpbXVsYXRlZF9oZWlnaHRzLCB+UEkoLngpWzFdKSwKICB1cHBlcl9waSA9IG1hcF9kYmwoc2ltdWxhdGVkX2hlaWdodHMsIH5QSSgueClbMl0pCikKYGBgCgojIDRIMgpgYGB7cn0KZGF0YSgiSG93ZWxsMSIpCmQgPC0gSG93ZWxsMQpkICU+JSAKICBmaWx0ZXIoYWdlIDwgMTgpICU+JSAKICBnZ3Bsb3QoYWVzKHggPSB3ZWlnaHQsIHkgPSBoZWlnaHQpKSArCiAgZ2VvbV9wb2ludCgpICsgZ2VvbV9zbW9vdGgobWV0aG9kPWxtKSArCiAgZ2VvbV9zbW9vdGgoY29sb3VyID0gIm9yYW5nZSIpCmBgYAoKYGBge3J9CmQgJT4lIAogIGZpbHRlcihhZ2UgPCAxOCkgJT4lIAogIG1hcChhbGlzdCgKICAgIGhlaWdodCB+IGRub3JtKG11LCBzaWdtYSksCiAgICBtdSA8LSBhICsgYiAqIHdlaWdodCwKICAgIGEgfiBkbm9ybShtZWFuID0gMTAwLCBzZCA9IDEwMCksCiAgICBiIH4gZG5vcm0obWVhbiA9IDAsIHNkID0gMTApLAogICAgc2lnbWEgfiBkdW5pZihtaW4gPSAwLCBtYXggPSA1MCkKICApLAogIGRhdGEgPSAuKSAtPgptb2RlbAoKcHJlY2lzKG1vZGVsKQpgYGAKCkEgb25lIHVuaXQgaW5jcmVhc2UgaW4gd2VpZ2h0IGluY3JlYXNlcyBoZWlnaHQgYnkgMi43MiwgdGhlcmVmb3JlIGEgMTAgdW5pdCBpbmNyZWFzZSB3aWxsIG1vdmUgaGVpZ2h0IGJ5IDI3LjIuCgpgYGB7cn0KIyBTYW1wbGUgdGhlIHBvc2V0ZXJpb3IKcG9zdGVyaW9yX3NhbXBsZXMgPC0gZXh0cmFjdC5zYW1wbGVzKG1vZGVsKQoKIyBDcmVhdGUgeC1heGlzCm5ld193ZWlnaHRzIDwtIHNlcShmcm9tID0gNCwgdG8gPSA0NSwgbGVuZ3RoLm91dCA9IDEwMCkKCiMgQ2FsY3VsYXRlIHRoZSBtZWFuCmZuX21lYW4gPC0gZnVuY3Rpb24od2VpZ2h0KSB7cG9zdGVyaW9yX3NhbXBsZXMkYSArIHBvc3Rlcmlvcl9zYW1wbGVzJGIgKiB3ZWlnaHR9CmZpdF9tZWFucyA8LSBwdXJycjo6bWFwKG5ld193ZWlnaHRzLCBmbl9tZWFuKQoKIyBDYWxjdWxhdGUgaGVpZ2h0cwpmbl9oZWlnaHQgPC0gZnVuY3Rpb24od2VpZ2h0KSB7CiAgcm5vcm0oCiAgICBuID0gbnJvdyhwb3N0ZXJpb3Jfc2FtcGxlcyksCiAgICBtZWFuID0gcG9zdGVyaW9yX3NhbXBsZXMkYSArIHBvc3Rlcmlvcl9zYW1wbGVzJGIgKiB3ZWlnaHQsCiAgICBzZCA9IHBvc3Rlcmlvcl9zYW1wbGVzJHNpZ21hKQp9CmZpdF9oZWlnaHRzIDwtIHB1cnJyOjptYXAobmV3X3dlaWdodHMsIGZuX2hlaWdodCkKCmRhdGFfZnJhbWUoCiAgd2VpZ2h0cyA9IG5ld193ZWlnaHRzLAogIG1lYW5zID0gbWFwX2RibChmaXRfbWVhbnMsIG1lYW4pLAogIG1lYW5zX2xvdyA9IG1hcF9kYmwoZml0X21lYW5zLCB+IEhQREkoLngsIC44OSlbMV0pLAogIG1lYW5zX2hpZ2ggPSBtYXBfZGJsKGZpdF9tZWFucywgfiBIUERJKC54LCAuODkpWzJdKSwKICBoZWlnaHRzID0gbWFwX2RibChmaXRfaGVpZ2h0cywgbWVhbiksCiAgaGVpZ2h0c19sb3cgPSBtYXBfZGJsKGZpdF9oZWlnaHRzLCB+IEhQREkoLngsIC44OSlbMV0pLAogIGhlaWdodHNfaGlnaCA9IG1hcF9kYmwoZml0X2hlaWdodHMsIH4gSFBESSgueCwgLjg5KVsyXSksCikgJT4lIApnZ3Bsb3QoYWVzKHggPSB3ZWlnaHRzKSkgKwogIGdlb21fbGluZShhZXMoeSA9IG1lYW5zKSkgKwogIGdlb21fcmliYm9uKGFlcyh5bWluID0gbWVhbnNfbG93LCB5bWF4ID0gbWVhbnNfaGlnaCksIGFscGhhID0gLjQpICsKICBnZW9tX3JpYmJvbihhZXMoeW1pbiA9IGhlaWdodHNfbG93LCB5bWF4ID0gaGVpZ2h0c19oaWdoKSwgYWxwaGEgPSAuMikgKwogIGdlb21fcG9pbnQoaW5oZXJpdC5hZXMgPSBGQUxTRSwgZGF0YSA9IGZpbHRlcihkLCBhZ2UgPCAxOCksCiAgICAgICAgICAgICBhZXMoeD13ZWlnaHQsIHkgPSBoZWlnaHQpLCBzaGFwZSA9IDEpICsKICB5bGFiKCJoZWlnaHQiKSArIHhsYWIoIndlaWdodCIpICsgdGhlbWVfYncoKSArCiAgdGhlbWUocGFuZWwuZ3JpZC5tYWpvciA9IGVsZW1lbnRfYmxhbmsoKSwgcGFuZWwuZ3JpZC5taW5vciA9IGVsZW1lbnRfYmxhbmsoKSkKYGBgCgojIEJSTVMKCmBgYHtyfQpkYXRhKEhvd2VsbDEpCmQgPC0gSG93ZWxsMQpkMiA8LSBkICU+JSAKICBmaWx0ZXIoYWdlID49IDE4KQpkMiAlPiUgCiAgc2tpbQpgYGAKCgpgYGB7cn0KbGlicmFyeShicm1zKQoKYjQuMSA8LQogIGJybShkYXRhID0gZDIsIGZhbWlseSA9IGdhdXNzaWFuLAogICAgICBoZWlnaHQgfiAxLAogICAgICBwcmlvciA9IGMocHJpb3Iobm9ybWFsKDE3OCwgMjApLCBjbGFzcyA9IEludGVyY2VwdCksCiAgICAgICAgICAgICAgICBwcmlvcihjYXVjaHkoMCwgMSksIGNsYXNzID0gc2lnbWEpKSwKICAgICAgaXRlciA9IDMxMDAwLCB3YXJtdXAgPSAzMDAwMCwgY2hhaW5zID0gNCwgY29yZXMgPSA0KQpwbG90KGI0LjEpCmBgYAoKRXh0cmFjdCB2YXJpYW5jZSBjb3ZhcmlhbmNlIG1hdHJpeC4KYGBge3J9CiMgRXh0cmFjdCBzYW1wbGVzIGZyb20gdGhlIHBvc3Rlcmlvcgpwb3N0IDwtIHBvc3Rlcmlvcl9zYW1wbGVzKGI0LjEpCiMgQ292YXJpYW5jZSBtYXRyaXgKY292KHBvc3RbLCAxOjJdKQpgYGAKCkxldHMgZXh0cmFjdCB0aGUgSERQSQpgYGB7cn0KbGlicmFyeSh0aWR5YmF5ZXMpCgpwb3N0ICU+JSAKICBzZWxlY3QoLWxwX18pICU+JSAKICBnYXRoZXIocGFyYW1ldGVyKSAlPiUgCiAgZ3JvdXBfYnkocGFyYW1ldGVyKSAlPiUgCiAgbWVkaWFuX2hkaSh2YWx1ZSkKYGBgCgpOb3cgbGV0cyBhZGQgYSBwcmVkaWN0b3IuCmBgYHtyIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9CmQzIDwtIGQyICU+JSAKICBtdXRhdGUod2VpZ2h0ID0gKHdlaWdodCAtIG1lYW4od2VpZ2h0KSkvc2Qod2VpZ2h0KSkgJT4lIAogIHNlbGVjdChoZWlnaHQsIHdlaWdodCkKCmI0LjMgPC0gCiAgYnJtKGRhdGEgPSBkMywgZmFtaWx5ID0gZ2F1c3NpYW4sCiAgICAgIGhlaWdodCB+IDEgKyB3ZWlnaHQsCiAgICAgIHByaW9yID0gYyhwcmlvcihub3JtYWwoMTU2LCAxMDApLCBjbGFzcyA9IEludGVyY2VwdCksCiAgICAgICAgICAgICAgICBwcmlvcihub3JtYWwoMCwgMTApLCBjbGFzcyA9IGIpLAogICAgICAgICAgICAgICAgcHJpb3IodW5pZm9ybSgwLCA1MCksIGNsYXNzID0gc2lnbWEpKSwKICAgICAgaXRlciA9IDQ2MDAwLCB3YXJtdXAgPSA0NTAwMCwgY2hhaW5zID0gNCwgY29yZXMgPSA0LAogICAgICBjb250cm9sID0gbGlzdChhZGFwdF9kZWx0YSA9IDAuOCwgCiAgICAgICAgICAgICAgICAgICAgIG1heF90cmVlZGVwdGggPSAxMCkpCgpwbG90KGI0LjMpCmBgYAoKRXhhbWluZSBjb3ZhcmlhbmNlCmBgYHtyfQpwYWlycyhiNC4zKQpgYGAKCkNyZWF0aW5nIHByZWRpY3Rpb25zCmBgYHtyfQpwb3N0IDwtIHBvc3Rlcmlvcl9zYW1wbGVzKGI0LjMpCgptdV9hdF81MCA8LQogIHBvc3QgJT4lIAogIHRyYW5zbXV0ZShtdV9hdF81MCA9IGJfSW50ZXJjZXB0ICsgYl93ZWlnaHQgKiA1MCkKCm11X2F0XzUwICU+JQogIGdncGxvdChhZXMoeCA9IG11X2F0XzUwKSkgKwogIGdlb21fZGVuc2l0eShzaXplID0gMCwgZmlsbCA9ICJncmV5NzUiKSArCiAgc3RhdF9wb2ludGludGVydmFsaChhZXMoeSA9IDApLCAKICAgICAgICAgICAgICAgICAgICAgIHBvaW50X2ludGVydmFsID0gbW9kZV9oZGksIC53aWR0aCA9IC45NSkgKwogIHNjYWxlX3lfY29udGludW91cyhOVUxMLCBicmVha3MgPSBOVUxMKSArCiAgbGFicyh4ID0gZXhwcmVzc2lvbihtdVsiaGVpZ2h0IHwgd2VpZ2h0ID0gNTAiXSkpICsKICB0aGVtZV9jbGFzc2ljKCkKCmBgYAoKTm93IHdlIGV4dHJhY3QgY29tcGxldGUgcHJlZGljdGlvbnMgaW5jbHVkaW5nIHRoZSBtZWFuIGFuZCBhZGRlZCB2YXJpYW5jZS4KYGBge3J9CndlaWdodF9zZXEgPC0gdGliYmxlKHdlaWdodCA9IHNlcShmcm9tID0gMjUsIHRvID0gNzAsIGJ5ID0gMSkpCgpwcmVkX2hlaWdodCA8LQogIHByZWRpY3QoYjQuMywKICAgICAgICAgIG5ld2RhdGEgPSB3ZWlnaHRfc2VxKSAlPiUKICBhc190aWJibGUoKSAlPiUKICBiaW5kX2NvbHMod2VpZ2h0X3NlcSkKICAKcHJlZF9oZWlnaHQgJT4lCiAgc2xpY2UoMTo2KQpgYGAKCmBgYHtyfQpnZ3Bsb3QocHJlZF9oZWlnaHQsIGFlcyh3ZWlnaHQsIEVzdGltYXRlKSkgKwogIGdlb21fcmliYm9uKGFlcyh5bWluID0gUTIuNSwgeW1heCA9IFE5Ny41KSwgZmlsbCA9ICJncmV5NzUiKSArCiAgZ2VvbV9saW5lKCkgKwogIHRoZW1lKHBhbmVsLmdyaWQgPSBlbGVtZW50X2JsYW5rKCkpCmBgYAoKV2UgY2FuIG92ZXJsYXkgdGhpcyB3aXRoIHRoZSBtZWFuIGludGVydmFsIHRvbwpgYGB7cn0KbXVfc3VtbWFyeSA8LQogIGZpdHRlZChiNC4zLCAKICAgICAgICAgbmV3ZGF0YSA9IHdlaWdodF9zZXEpICU+JQogIGFzX3RpYmJsZSgpICU+JQogIGJpbmRfY29scyh3ZWlnaHRfc2VxKQoKZ2dwbG90KG11X3N1bW1hcnksIGFlcyh3ZWlnaHQsIEVzdGltYXRlKSkgKwogIGdlb21fcmliYm9uKGFlcyh5bWluID0gUTIuNSwgeW1heCA9IFE5Ny41KSwgZmlsbCA9ICJncmV5NzUiKSArCiAgZ2VvbV9saW5lKCkgKwogIHRoZW1lKHBhbmVsLmdyaWQgPSBlbGVtZW50X2JsYW5rKCkpCmBgYAoKYGBge3J9CmlubmVyX2pvaW4ocHJlZF9oZWlnaHQsIG11X3N1bW1hcnksIGJ5ID0gIndlaWdodCIsIHN1ZmZpeD1jKCJfaGVpZ2h0IiwgIl9tZWFuIikpICU+JSBzdHIKYGBgCgpgYGB7cn0KaW5uZXJfam9pbihwcmVkX2hlaWdodCwgbXVfc3VtbWFyeSwgYnkgPSAid2VpZ2h0Iiwgc3VmZml4PWMoIl9oZWlnaHQiLCAiX21lYW4iKSkgJT4lIAogIGdncGxvdChhZXMoeCA9IHdlaWdodCwgeSA9IEVzdGltYXRlX2hlaWdodCkpICsKICBnZW9tX3JpYmJvbihhZXMoeW1pbiA9IFEyLjVfaGVpZ2h0LCB5bWF4ID0gUTk3LjVfaGVpZ2h0KSwgZmlsbCA9ICJncmV5MjUiKSArCiAgZ2VvbV9yaWJib24oYWVzKHltaW4gPSBRMi41X21lYW4sIHltYXggPSBROTcuNV9tZWFuKSwgZmlsbCA9ICJncmV5NDUiKSArCiAgZ2VvbV9saW5lKCkgKwogIHRoZW1lKHBhbmVsLmdyaWQgPSBlbGVtZW50X2JsYW5rKCkpCmBgYAoKCgoKCgoKCgoKCg==

Chapter 4 - Linear Models